## --------------------------------------- ##
## simulation setup: 
##  delta = n / p
##  rho = k / n
## --------------------------------------- ##

## --------------------------------------- ##
## cluster setup 
## --------------------------------------- ##
## i <- as.numeric(Sys.getenv("SLURM_ARRAY_TASK_ID")) # this should be 1--50
for(i in 1:length.out){

  cl <- makeCluster(detectCores()/2)
  registerDoParallel(cl)
  

## draw error term and n 
set.seed(38193)
n    <- floor(p * delta[i])
eta  <- rnorm(n)
etaz <- sigma * eta 

# doRNGseed(1234)
sim_out <- foreach(j = 1:length(rho), .combine = "rbind", .packages= c("BridgeChange", "glmnet", 'mvnfast')) %dorng% ({
    
    cat("Now at i = ", i, " j = ", j ,"\n")
    k <- ceiling(n * rho[j])

    l2_error <- save_bridge <- save_ridge <- save_en  <- rep(NA, n_iter)
    for (iter in 1:n_iter) {
      rm(y); rm(X);
        ## gen data ***************************************************************
        beta <- c(rnorm(k, 0, 3), rep(0, p - k))
        Xfull <- as.matrix(faux::rnorm_multi(n = n,  mu = rep(0, p), r = 0.7))
        ## rmvn(n, rep(0, p), Sigma)
        beta <- c(rnorm(k, 0, 3), rep(0, p - k))
        yfull <- Xfull %*% beta + etaz

        n_out <- floor(n * n_miss)
        id_use <- sample(c(rep(0, n_out), rep(1, n - n_out)))
        y <<- yfull[id_use == 1]; X <<- Xfull[id_use == 1, ]

        ## fit lasso ***************************************************************
        # out <- rlasso(y = y, x = X, intercept = FALSE)
        # l2_error[iter] <- sqrt(sum((coef(out) - beta)^2)) / sqrt(sum(beta^2))
        out <- cv.glmnet(y = y, x= X, alpha = 1, nfolds = 3)
        est_error <- yfull[id_use ==0] - Xfull[id_use == 0,] %*% as.vector(coef(out, s = "lambda.min")[-1,])
        l2_error[iter] <- mean(est_error^2)

        ## fit SparseChange ********************************************************
        out1 <- BridgeChangeReg(y~X, mcmc= 100, burn = 100, thin=1, verbose=0, n.break = 0)
        beta.bridge <- coef_bridge(out1)
        save_bridge[iter]  <- mean((yfull[id_use == 0] - Xfull[id_use == 0,] %*% beta.bridge)^2)

        ## fit ridge ********************************************************
        rg_out <- cv.glmnet(y = y, x= X, alpha = 0, nfolds = 3)
        est_error <- yfull[id_use ==0] - Xfull[id_use == 0,] %*% as.vector(coef(rg_out, s = "lambda.min")[-1,])
        save_ridge[iter] <- mean(est_error^2)

        ## fit elastic net **************************************************
        en_out <- cv.glmnet(y = y, x= X, alpha = 0.5, nfolds = 3)
        est_error <- yfull[id_use ==0] - Xfull[id_use == 0,] %*% as.vector(coef(en_out, s = "lambda.min")[-1,])
        save_en[iter] <- mean(est_error^2)

    }

    mse.lasso   <- median(l2_error, na.rm = TRUE)
    mse.bridge  <- median(save_bridge, na.rm = TRUE)
    mse.ridge   <- median(save_ridge, na.rm = TRUE)
    mse.elastic <- median(save_en, na.rm = TRUE)


    out <- c(mse.lasso, mse.bridge, mse.ridge, mse.elastic)
    out
})

    cat("\nCurrent simulation is ", i, " iteration.\n")

    saveRDS(sim_out, file = paste("./nochange/corr07/res_cv/sim_cv_mse_corr07_", i, ".rds", sep = ''))
    
    
    stopCluster(cl)
    stopImplicitCluster()
}
